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Abstract. We consider a system of N individuals consisting of S species that interact 
pairwise: x m + xg —> 2x m with arbitrary probabilities p m . With no spatial structure, 
QJ ' the master equation yields a simple set of rate equations in a mean field approximation, 

£h the focus of this note. Generalizing recent findings of cyclically competing threc- 

and four-species models, we cast these equations in an appealingly simple form. As 
a result, many general properties of such systems are readily discovered, e.g., the 
. major difference between even and odd S cases. Further, we find the criteria for the 

existence of (subspaces of) fixed points and collective variables which evolve trivially 
(exponentially or invariant). These apparently distinct aspects can be traced to the 
null space associated with the interaction matrix, p l m . Related to the left- and right- 
zero-eigenvectors, these appear to be "dual" facets of the dynamics. We also remark 
on how the standard Lotka-Volterra equations (which include birth/death terms) can 



Q , be regarded as a special limit of a pairwise interacting system. 

Introduction. Population dynamics is a venerable subject, dating back two centuries 
to Malthus, Verhulst, Lotka, Volterra, and many others [H [21 [3]. Nonetheless, new and 
interesting phenomena are continually being discovered. For example, many studies 
of cyclic competition between 3 species (with no spatial structure, e.g., a well- mixed 
system) attracted considerable recent attention [3]. In fact, in 2009, Science Daily 
popularized this topic[5] by branding it "Survival of the Weakest." We extended this 
investigation to a system with 4 species[6], which displayed no such counter-intuitive 
behavior. Instead, we found an intuitively understandable principle which underpins all 
systems with cyclically competing species, namely, "The prey of the prey of the weakest 
is least likely to survive." In the case of cyclically competing 3 species, the prey of one's 
prey is also one's predator. Thus, its demise is indeed "good news" for the weakest and 
leads to the eye-catching headline. In this short note, we considered a wider range of 
systems of S species interacting pairwise, with arbitrary rates. Focusing only on a mean 
field description (i.e., rate equations), we find remarkable general properties, such as 
fixed points, invariant manifolds, collective variables with simple time dependence, as 
well as a (possibly new) form of "duality." We begin by specifying the individual based 
stochastic model, from which our MF approximation is derived. This note is devoted 
only to the properties of solutions to the MF equations, however. 

Individual based model. Consider a system with N individuals, each being a member 
of one of S species. Let us denote the species by x m , with m = 1, S, and the number 
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of individuals of each by N m . We allow only pairwise interactions, i.e., 

x m + xt 2x m . (1) 

where p m are arbitrary probabilities for a "predator" x m to consume a "prey" xe- Note 
that, if we wish to model bi-directional interactions, then these p's represent the net 
consumption of the dominant species. Thus, there are at most S (S — 1) /2 such positive 
quantities. When two individuals encounter, the role of each is well defined: one is the 
predator, the other is the prey. To emphasize, we illustrate with p\ = 0.7, and both 
X\ + x 5 and x 5 + X\ becomes 2x\ with probability 0.7 (and unchanged with probability 
0.3). Defining an update of our system as randomly choosing a pair and letting them 
interact, we see that the iV m 's change by ±1 or with N = Yl m remaining a constant 
for all times. Depending on the p^'s, there are at least S absorbing states. 

The appropriate quantity for describing this stochastic evolution is, of course, 
P {{N m } ; t), the probability to find the system with the set {N m } after updating it r 
times from some given initial P ({N m } ; 0). The change in one step, P ({N m } ; r + 1) — 
P ({N m } ; r), is given by 2/N (N — 1) times 

Y,P e n l(N n ~ 1) (N t + 1)P ({N m - 5 mn + 5 mi } ; r) - N n N e P ({N m } ; r)] (2) 

(n,e) 

where the sum is over only predator-prey pairs. Finally, to include the standard form 
of a Lotka-Volterra model, we should add birth/death probabilities 

x m 2x m , x e (3) 

as well. Except for the last section, in which we will comment on this addition, we 
consider only pair interactions in the remainder of this note. 

Rate equations and general properties of their solutions. To proceed, we exploit 
the standard MF approximation. Taking the large N limit, rescaling time by t = r/N 
(which becomes a continuous variable), and replacing averages 

o (r) = (o) T (W) p ; r ) 

{N m } 

of products of N m by the product of the averages, we arrive at the rate equations for 
the fractions 



namely, 



Here, 



(X m ) T = (N m )jN (4) 



dX m = = X m ^2 k m eX E . (5) 



Km = -hm = 2p l ri 
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can be regarded as the elements of an antisymmetric matrix K. Although eqn. (jSJ) is 
nonlinear, its special form allows us to write a quasi-linear equation 0. Using vector 



notation, it is 



d\U) = K\X) (6) 

where the elements of \U) are 

U m = \nX m (7) 

The rest of this note is devoted to the consequences of eqns. (jSJ) or (jSJ), apart from the 
obvious, strict conservation law: Yli Xg = 1. Denoting a row vector with all its elements 
being unity by (1|, this constraint reads 

(1|X) = 1. (8) 

For physical systems, we also demand Xi > 0, so that the evolution of interest here 
takes place in an S — 1 simplex, i.e., a hypertetrahedron in S — 1 dimensions. 

If A is an eigenvalue of K, so is —A, since K is antisymmetric. Thus, if S is odd, 
K must have at least one zero eigenvalue. In case there is only one, let us denote its 
associated left and right eigenvectors by: 

(C|K = K|C) = 0. (9) 

Note that, in general, a left eigenvector is not simply related to the transpose of the right 
eigenvector. However, since — K|C) = = (— K|C)) T = (|£)) T IK, (£| is the transpose of 
These two eigenvectors play "dual" roles in the following sense. Obviously, \() is a 
fixed point of eqn. (J2]) and, if every ( m is non-negative %, then 

X* m = Cm^J2^ ( 10 ) 

represents a fixed population [JJ. If some are negative, then the system has no fixed 
populations (except for absorbing states). Meanwhile, we also have d (£| U) = 0, so 
that 

n = \[{x m f m (ii) 

m 

is a constant in the evolution. In the three cyclically competing species case [I], the 
invariant R = A kb B kc C ka can be readily related to 1Z. When all £ m 's are positive, we 
see that these invariant manifolds (i.e., 7£'s given by various initial values of X m ) are 
closed S — 2 dimensional subspaces within our S — 1 simplex - e.g., the closed loops in 

| A special system, with a form similar to (|5|), can be solved exactly [2]. Instead of a matrix, k m £, 
there is just a constant vector gg. Such a system is suitable for species competiting for the same pool 
of resources, which gets depleted via the combination ^2 f gpXi. I thank Zoltan Toroczkai for pointing 
this case out to me. 

§ More precisely, we need CnCm to be non-negative for all n, m. 

|| Strictly, fixed points with some components being zero are also uninteresting, since an extinct speices 
will remain so for ever. For simplicity, readers may restrict their attention to points within the S — 1 
simplex. 
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[I]. For S > 3, we can be certain of the system evolving along orbits that lie in such 
a manifold, but we have not shown whether these orbits are open or closed. In either 
case, the long time average of X m (t), defined by 



1 f T 

X m = lim — / X m (t) dt 



is precisely X^. The proof is simple: Integrate eqn. (jBJ), recognize the finiteness of 
\nX m , and invoke the uniqueness of \(). Note that 1Z does change in a stochastic 
model, as these orbits and the fixed point are neutrally stable [7]. On the other hand, 
if some £ m 's are negative, 1Z is still meaningful and typically determine how pairs of 
species must vanish in the long time limit. A good example is the three species case 
with say, A preying on C with rate k c (instead of being consumed by C). Then, C must 
decrease monotonically. The (related) invariant is R = A kb B~ kc C ka , so that B and C 
must vanish together, with C — > \RB c \ a , leaving A as the sole survivor. Finally, 
if the spectrum of IK contains v > 1 zeros, then there will be v invariants and, "dual" 
to those, a v — 1 dimension subspace of fixed points. Many of the conclusions can be 
readily generalized. It would be interesting to explore further details of these physical 
subspaces, especially if some of these eigenvectors have negative components. 

The case of even S displays even richer behavior. First, it is possible that none of 
the eigenvalues of K vanishes. In that case, it has a unique inverse, so that (EJ can be 
written as 

OK' 1 \U) = \X) (12) 
so that d (1| K- 1 \U) = 1. Defining 

Vm = J2(K-i) em (13) 

i 

we see that the collective variable 

Q = \[(X m T™ (14) 

m 

evolves trivially: 

Q(t) = Q(0)e t (15) 

At first glance, this mathematical expression appears very "rigid," as it predicts an 
exponentially increasing Q under all cirmstances. With a little reflection, a more 
"flexible" physical interpretation emerges: Whether a species wins or loses is associated 
with the signs of the rfs. Specifically, at least one of the species associated with negative 
77's must vanish at late times. 

Let us turn to some simple examples to illustrate these findings. In the most trivial 
case of S = 2, let x\ consume £2 with rate k > 0. Thus, X2 must vanish at large 
t, consistent with Q being (X 1 /X 2 ) 1 ^ fc . Indeed, the full solution is given simply by 
X 1 /X 2 = e kt Xi/X 2 \ t _ . Since X\ + X 2 — 1, the explicit solution follows from algebra. 
A particularly appealing form exploits the symmetry X x X 2 © t <^ — t and displays 
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the symmetric sigmoid associated with the evolution from X 2 (—00) = 1 to X% (00) = 1. 
Defining 

kt 

a (t) = tanh — (16) 

we easily find Xi(t) — X 2 (t) = o-(t — to), with i given by the initial condition: 
(2/fc) tanh -1 [X 2 (0) — Xi (0)]. Thus, the explicit solution is 

^.,w = 1±g( 2 t " fo) . (17) 

At the next level, S = 4, we have the special case of four cyclically competing 
species where the variable Q = A kb+kc B~ kc ~ kd C kd+ka D~ ka ~ kb can be readily related 
to Q. Further, this system is an excellent example of a large class of even S systems, 
namely, ones with two "opposing teams." Specifically, suppose each "team" consists 
of s = S/2 species - {yj} ,{zj} ; j = 1, s - with the understanding that species on 
the same team do not interact. Then, all the interactions can be specified by the s 2 
probabilities associated with encounters between y^ and Zj. Note that the outcome can 
be either 2y^ or 2zj, depending on which player is stronger. Note also that the absorbing 
states of the stochastic model form two s— 1 simplices. To proceed, we let (Yl, Zj) denote 
the fractions of the species (yj, zf) and arrange the column vector according to 

Then, K takes the form 





itt 







(19) 



Here M is a s x s matrix, the elements of which are associated with the rate of yi 
consuming (or being consumed by) Zj. The antisymmetry of K is "built in" and, unlike 
its spectrum, there are no constraints on the spectrum of M. Focusing first on non- 
singular M's, we define 

A = det M ^ 

(generalization of A in [B]) and 

W = AM" 1 (20) 

the elements of which are the appropriate co-factors of M. These considerations 
allows us to find a non-trivial invariant manifold - a straight line joining the point 
W T |1) / (1| W T |1) in the space of \Y) fij with the point W |1) / (1| W |1) in the space of 
\Z). Indeed, the evolution on this line is precisely the same as the S = 2s = 2 case, 
dll, i.e., 

l/(l + ,)W-|l> 
17 2 \ (1 -<t)W|1) 1 



% Here, |1) = . We also assume W |1) and W |1) have positive elements. The cases where this is 
not true deserve further study. 
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provided k in (TTSj) is replaced by A and an appropriate t is taken. The proof is 
straightforward, since d In [7/ (t)] =dlnf (t) for any constant 7, so that 

Meanwhile, we have 
K\X) = 

But <91n(l±(x) = ±do-/(l ±0-) = ±A(1=fo-)/2, so that d\U) = K\X) is indeed 
satisfied with a = tanh(Ai/2). The result (12~T1) is appealing, since A = det M is 
precisely k in case s = 1. In [9], the case of s = 2 will be presented in explicit detail. 
Note that we do not assume a sign for A here. If A > 0, team {yj} will "win," i.e., the 
system ending at \Z) = (and \Y) oc W T |1)). Similarly, if A < 0, the system will end at 
\Z) oc W Apart from this special line, we can draw similar conclusions based simply 
on the sign of A, namely, which team "wins" and which "loses." We should caution the 
reader that this way of displaying extinction is restricted to systems with two teams, 
while eqns. (114)15p are valid in general. 

Next, we turn to systems with singular K's, i.e., neutrally stable ones. The null 
space of IK must be even dimensional, a number we denote by 2/i (integer As in the 
odd S case, the associated left- and right-eigenvectors play "dual" roles in that they 
can be used to construct 2/i invariants as well as a 2/i — 1 dimensional subspace of fixed 
points (which may or may not be physical). A good illustration is the special case of 
S = 2s cyclically competing species, for which a straightforward (though tedious) proof 
leads to /i being unity always [9]. More generally, for systems with two opposing teams, 
we can exploit (119"]) and show explicit expressions for the subspace of fixed points and 
invariants. Clearly, the methods for odd S applies here: Define the zero eigenvectors of 
M by 

(Ca|M = M|C«) = 0. (22) 
with a — 1, jJL. Note that, unlike K, M is not antisymmetric (or symmetric) in general, 



Ca) (i-e., the tranpose of (C 



so that M T is not simply related to M. Thus, 
are linearly independent, typically. The subspace of fixed points is spanned by the 2\i 
vectors 




and \Ca) 



(23) 

In case all the elements of these vectors are non-negative, this subspace will have a 
non-trivial intersection with the physical region: X m > and ^2 m X m = 1. Then, the 
subspace of physical fixed points is 2/i — 1 dimensional. In all cases, the 2/i invariants 
defining orbits of neutral stability can be written: 

s s 
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showing again the "dual" role played by the zero-eigenvectors. 

So far, we considered systems with non-singular and singular K's separately. A 
unified approach can be formulated, taking advantage of the cofactors, W, being well 
defined even if A = and M -1 does not exist. Starting with eqns. ( 16|7|) . we multiply by 








(25) 



The right side is simply A \X) by virtue of eqn. ( 1201) . whether A = or not. Taking 
the inner product with (1|, the right side reduces to A. Turning to the left side, let us 
denote the elements of (1| W and (1| W T by £j and Xj respectively. These allow us to 
define two other collective variables 

/••-n^r ; G^niz^ (26) 

3 3 

From these, we form 

Q(t) = ^ = Q (0) e At . (27) 

which is readily seen as Q but is more "flexible" - in that Q can increase or decrease 
with t. Similarly, it is straightforward to show that, as A — > 0, F and G reduces to 1Z 
and TZ in ( I24jl . In this sense, the variables ( 1261271) are superior to ( l24lfT4l) . 

Concluding remarks. In this brief note, we report findings of general properties 
associated with the rates equations for a system of S species that compete pairwise. In 
some cases (all odd S and some even S), subspaces of fixed points appear to be "dual" to 
invariant manifolds. In the absence of this behavior, we find a collective variable which 
evolves exponentially. We conclude with a few remarks on possible future research, 
apart from the issues already raised above. 

In a standard Lotka-Volterra model of predator-prey interactions, individuals can 
be born and can die without the interaction with those of another species, as in (j3J). 
Thus, an equation like (jSJ) should contain linear terms 

dN m = lm N m + N m J2 k meN e (28) 

e 

where 7 m is b m or — d m . Of course, we can re-analyze our problem with these additions. 
However, it is possible to regard this problem as a special limit of an S + 1 species 
system. First, note that 7 m are generally chosen to be similar to k m £, so that the typical 
iV m 's reflect the numbers in nature (thousand, million, etc.). For example, in the S = 2 
case, the typical levels of the predator and prey are d/k and b/k, respectively. Next, 
introduce N individuals of species x , let them interact with x m ^ according to 

x m + x 2x m ; x m + x 2x (29) 

and consider the limit of N,N — > oo. Thus, X = O (1) and X m ^ = O (1/N). Now, 
consider a generalized eqn. (J5]): 

dX m = X m L m iXe (30) 
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with m, £ = 0, 1, S and 

L ml = Nk mt m,£^0 

Keeping lowest order terms in each equation, we find a consistent limit for (J28l) . i.e., 

<9V = + 0(1/ N) (31) 
(which admits X = 1 + 0(1/N)) and, for m^O, 

«9V m = 7m X m + X m kmeNi + 0{l/N 2 ) (32) 

i 

In other words, the evolution of an ordinary Lotka-Volterra system takes place in a tiny 
"corner" of the enlarged configuration space. In particular, the familiar fixed point and 
closed orbits in the 2 species Lotka-Volterra system (d/k,b/k) are just the ones in the 
case of 3 cyclically competing species. 

Since all realistic evolution is stochastic, this study should be extended, along the 
lines in [7], for example. The invariants here are clearly related to slow variables and 
identifying the fast variables will be helpful. For systems with two teams, expression 
(fl9l) behooves us to explore a possible simplectic structure in our dynamics, even though 
we must deal with difficulties associated with non-linearities. In several special cases, 
the implications of Nambu dynamics for this kind of evolution was studied [H]. Since 
that approach requires S — 1 invariants ( "Hamiltonians" ) , we can easily see if it can be 
exploited for any given set {p^}, by checking whether the spectrum of K has S — 2 
zero's. Beyond such possibilities, an entirely new vista of systems with competing species 
awaits us, e.g., ones with non-trivial spatial structures, ones on complex networks, and 
ones with inhomogeneous rates. Clearly, to reflect the enormous range of phenomena in 
nature, we can create many more models and can expect many more interesting results. 
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